disp('Plotting a few bonus figures...')
disp(sprintf('\n In file ''%s''',mfilename))
disp('These are some rough figures...file could use a lot of editing...')

warning('off','MATLAB:legend:IgnoringExtraEntries')



ConstraintOptions= {'dFdPwaterchange','rFrac_crit','Water','Tmp','U','Bulk Comp','FirstMeltPb','MeltExtractedPb','LastMeltPb','Gar2SpTransition','Sp2PlagTransition','aluminumOutPb','aluminumOutFrac','cpxOutPb','cpxOutFrac','Tdiff','quick crustal thickness','LastAdiabaticPressure','LastAdiabaticPressure2','LastAdiabaticPressure3','BrittleDuctile600Pres','ct2_ininc(1)','percentGar','percentSpin','percentPlag','percentGar_onaxis','percentSpin_onaxis','percentPlag_onaxis','FP_Full','FP_OnAxis','Fmax_OnAxis','Fmax_OnAxis_NotExp','230Th ZI FULL','226Ra ZI FULL','230Th CI FULL','226Ra CI FULL','maxNF 230Th ZI','maxNF 226Ra ZI','maxNF 230Th CI','maxNF 226Ra CI','TYPE','percentPool','distance_ininc','Variable Pool Fmean','Vpool %gar','Vpool  %sp','Vpool %plag','gar km','sp km','plag km','totalkm','Pb','T','Ts','Ln_cumulative','Lnexp_cumulative','StabilityField ','FN or FP','ZetaGar','ZetaSpinel','ZetaPlag','ColumnNum','variablypooledLnExp','variablypooledPb','variablypooledTs'};

[a,xxxvalue]=ismember('MeltExtractedPb',ConstraintOptions);
MeltExtractedPb = MASTER_SM_INFO(:,xxxvalue);

[a,xxxvalue]=ismember('FirstMeltPb',ConstraintOptions);
FirstMeltPb  = MASTER_SM_INFO(:,xxxvalue);

[a,xxxvalue]=ismember('Water',ConstraintOptions);
initialwatercontent  = MASTER_SM_INFO(:,xxxvalue);

[a,xxxvalue]=ismember('percentGar',ConstraintOptions);
Garnet_proportion  = MASTER_SM_INFO(:,xxxvalue);

[a,xxxvalue]=ismember('LastMeltPb',ConstraintOptions);
LastMeltPb  = MASTER_SM_INFO(:,xxxvalue);

markershape='ok';
MarkerColorHere='k';




%%
figure()
hold on
subaxis(2,2,1)
axis square
hold on
plot(initialwatercontent, FirstMeltPb,markershape,'MarkerSize',12,'MarkerFaceColor',MarkerColorHere)

subaxis(2,2,2)
axis square
hold on
plot(initialwatercontent, MeltExtractedPb,markershape,'MarkerSize',12,'MarkerFaceColor',MarkerColorHere)

subaxis(2,2,3)
axis square
hold on
plot(initialwatercontent, Garnet_proportion,markershape,'MarkerSize',12,'MarkerFaceColor',MarkerColorHere)

subaxis(2,2,4)
axis square
hold on
plot(initialwatercontent, LastMeltPb,markershape,'MarkerSize',12,'MarkerFaceColor',MarkerColorHere)



subaxis(2,2,1)
hold on
ylabel('First Melt Pressure (kbar)')
xlabel('Water in Source (ppm)')
axis([0 1000 10 60])
box on
ticklengthgood = 3*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
set(gca, 'Xtick',[0:100:1000])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 50;
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
title = 'Water_v_FirstMeltPressure';
set(gcf,'name',title)
set(gca,'YDir','Reverse')


subaxis(2,2,2)
ylabel('First Extracted Melt Pressure (kbar)')
xlabel('Water in Source (ppm)')
axis([0 1000 10 60])
axis square
box on
ticklengthgood = 3*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
set(gca, 'Xtick',[0:100:1000])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 50;
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
title = 'Water_v_FirstExtractedMeltPressure';
set(gcf,'name',title)
set(gca,'YDir','Reverse')

subaxis(2,2,3)
axis([0 1000 0 .25])
axis square
ylabel('Fraction Garnet in Pooled Melt')
xlabel('Water in Source (ppm)')
set(gca, 'Xtick',[0:100:1000])
box on
ticklengthgood = 3*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 50;
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
%increment = .005;
increment = .01;
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
title = 'Water_v_WGP'; set(gcf,'name',title)


subaxis(2,2,4)
ylabel('Last Melt Pressure (kbar)')
xlabel('Water in Source (ppm)')
axis([0 1000 0 30])
box on
ticklengthgood = 3*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
set(gca, 'Xtick',[0:100:1000])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 100;
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
title = 'Water_v_LastMeltPressure';
set(gcf,'name','Water in Source v Things')
set(gca,'YDir','Reverse')
set(gcf,'position',[258 1 837 704])



%%
figure()
hold on
set(gcf,'Position',[179 373 1162 531])
subaxis(1,2,1,'Margin',.03,'Spacing',.02,'PaddingBottom',.08)
hold on
axis square
ColorsHERE
zeta_N = MASTER_INSTANTS_INFO(:,39+4);
zeta_N(zeta_N==1)=1;
zeta_N(zeta_N~=1)=0;

%plot(MASTER_INSTANTS_INFO(:,19), MASTER_INSTANTS_DISEQ(:,1),'k.')

% plot(zeta_N, MASTER_INSTANTS_DISEQ(:,1),'k.')
% plot(MASTER_POOLS_INFO(:,19), MASTER_POOLS_DISEQ(:,1),'ko','MarkerSize',12)
TallOptions = [1300 1350 1400 1450];
Tcolors = {blue green yellow red};

for ss = [4 3 2 1]
    N_index = find(MASTER_INSTANTS_INFO(:,4)==TallOptions(ss));
    P_index = find(MASTER_POOLS_INFO(:,4)==TallOptions(ss));
    plot(zeta_N(N_index), MASTER_INSTANTS_DISEQ(N_index,1),'k.','color',Tcolors{ss})
    plot(MASTER_POOLS_INFO(P_index,19+4), MASTER_POOLS_DISEQ(P_index,1),'ko','MarkerSize',12,'color',Tcolors{ss})
end



legend('Zero Ingrowth Full NF','Zero Ingrowth Pool','Location','Best')
axis([-.1 1.1 .2 2])
set(gca,'fontsize', 20,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
%set(gca,'Ytick',[0:5:60])
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
ticklengthgood = 2*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
%set(gca,'YDir','Reverse')
xlabel('\zeta')
ylabel('(230Th/238U)')
%set(gcf,'Position',[360 84 346 614])
box on
set(gcf,'Name','(230Th/238U) v WGP')

subaxis(1,2,2)
hold on
axis square
% plot(zeta_N, MASTER_INSTANTS_DISEQ(:,3),'k.')
% %plot(MASTER_INSTANTS_INFO(:,19), MASTER_INSTANTS_DISEQ(:,3),'k.')
% plot(MASTER_POOLS_INFO(:,19), MASTER_POOLS_DISEQ(:,3),'ko','MarkerSize',12)

for ss = [4 3 2 1]
    N_index = find(MASTER_INSTANTS_INFO(:,4)==TallOptions(ss));
    P_index = find(MASTER_POOLS_INFO(:,4)==TallOptions(ss));
    plot(zeta_N(N_index), MASTER_INSTANTS_DISEQ(N_index,3),'k.','color',Tcolors{ss})
    plot(MASTER_POOLS_INFO(P_index,19+4), MASTER_POOLS_DISEQ(P_index,3),'ko','MarkerSize',12,'color',Tcolors{ss})
end


legend('Complete Ingrowth NF', 'Complete Ingrowth Full Pool','Location','Best')

axis([-.1 1.1 .2 2])
set(gca,'fontsize', 20,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
%set(gca,'Ytick',[0:5:60])
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
ticklengthgood = 2*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
%set(gca,'YDir','Reverse')
xlabel('\zeta')
ylabel('(230Th/238U)')
%set(gcf,'Position',[360 84 346 614])
box on
set(gcf,'Name','230Th excess v WGP')

%% NEW
figure()
hold on
set(gcf,'Position',[179 373 1162 531])
subaxis(1,2,1,'Margin',.03,'Spacing',.02,'PaddingBottom',.08)
hold on
axis square

zeta_N = MASTER_INSTANTS_INFO(:,37+4);
% zeta_N(zeta_N==1)=1;
% zeta_N(zeta_N~=1)=0;

%plot(MASTER_INSTANTS_INFO(:,19), MASTER_INSTANTS_DISEQ(:,1),'k.')

% plot(zeta_N, MASTER_INSTANTS_DISEQ(:,1),'k.')
% plot(MASTER_POOLS_INFO(:,19), MASTER_POOLS_DISEQ(:,1),'ko','MarkerSize',12)
% TallOptions = [1300 1350 1400 1450];
% Tcolors = {'b' 'g' 'y' 'r'};


for ss = [4 3 2 1]
    N_index = find(MASTER_INSTANTS_INFO(:,4)==TallOptions(ss));
    P_index = find(MASTER_POOLS_INFO(:,4)==TallOptions(ss));
    plot(100.*zeta_N(N_index), MASTER_INSTANTS_DISEQ(N_index,1),'k.-','color',Tcolors{ss})
    plot(100.*MASTER_POOLS_INFO(P_index,22+4), MASTER_POOLS_DISEQ(P_index,1),'ko','MarkerSize',12,'color',Tcolors{ss})
end



legend('Zero Ingrowth Full NF','Zero Ingrowth Pool','Location','Best')
axis([0 20 .2 2])
set(gca,'fontsize', 20,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
%set(gca,'Ytick',[0:5:60])
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 1; %findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
ticklengthgood = 2*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
xtickformat(gca, 'percentage');
%set(gca,'YDir','Reverse')
xlabel('F')
ylabel('(230Th/238U)')
%set(gcf,'Position',[360 84 346 614])
box on
set(gcf,'Name','(230Th/238U) v F')

subaxis(1,2,2)
hold on
axis square
% plot(zeta_N, MASTER_INSTANTS_DISEQ(:,3),'k.')
% %plot(MASTER_INSTANTS_INFO(:,19), MASTER_INSTANTS_DISEQ(:,3),'k.')
% plot(MASTER_POOLS_INFO(:,19), MASTER_POOLS_DISEQ(:,3),'ko','MarkerSize',12)

for ss = [4 3 2 1]
    N_index = find(MASTER_INSTANTS_INFO(:,4)==TallOptions(ss));
    P_index = find(MASTER_POOLS_INFO(:,4)==TallOptions(ss));
    plot(100.*zeta_N(N_index), MASTER_INSTANTS_DISEQ(N_index,3),'k.-','color',Tcolors{ss})
    plot(100.*MASTER_POOLS_INFO(P_index,22+4), MASTER_POOLS_DISEQ(P_index,3),'ko','MarkerSize',12,'color',Tcolors{ss})
end


legend('Complete Ingrowth NF', 'Complete Ingrowth Full Pool','Location','Best')

axis([0 20 .2 2])
set(gca,'fontsize', 20,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
%set(gca,'Ytick',[0:5:60])
ax = gca;
ax.XAxis.MinorTick = 'on';
increment = 1; %findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = .1; %findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
ticklengthgood = 2*[0.01 0.025];
set(gca,'ticklength',ticklengthgood)
xtickformat(gca, 'percentage');
%set(gca,'YDir','Reverse')
xlabel('F')
ylabel('(230Th/238U)')
%set(gcf,'Position',[360 84 346 614])
box on
set(gcf,'Name','230Th excess v F')

%%
% trace element plots

clear layerplotnum
figure('units','normalized','outerposition',[0.3935 0.3238 0.3923 0.5314])


hold on
box on
axis square
set(gca,'fontsize', 15,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
for layerplotnum = 1:runnumber
    %         if ischar(layerplot_color{layerplotnum})==1
    %             layerplot_color_here = char(layerplot_color(layerplotnum));
    %         else
    %             layerplot_color_here = layerplot_color{layerplotnum};
    %         end
    layerplot_color_here = 'k';
    layerplot_marker = '';
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)
    
    CPX=eval(sprintf('CS_minerals_OUT_%g(:,:,1)', layerplotnum));
    
    %CPX=eval(sprintf('CS_minerals_OUT_%g', layerplotnum));
    %CPX_1 = CS_minerals_OUT_1(:,:,1);
    % Ti is 4
    % Tb is 26
    plot(CPX(:,8), CPX(:,4),'.k')
    hold on
    plot(CPX(end,8), CPX(end,4),'*r')
    
    
    %plot(eval(sprintf('CumFrac_%g', layerplotnum)),eval(sprintf('test_%g(:,10)', layerplotnum)),char(layerplot_marker(layerplotnum)),'Color',layerplot_color_here,'linewidth',3)
end
xlabel('Zr')
ylabel('Ti')
set(gca,'ticklength',2*[.02 .05])
%axis([.01 100 100 10000])
set(gca,'yscale','log')
set(gca,'xscale','log')
set(gcf,'Name','ZrTi_in_CPX')


clear layerplotnum
figure('units','normalized','outerposition',[0.3935 0.3238 0.3923 0.5314])
%FigureTitle = sprintf('Ti Zr Tmp: %d %cC, %s', Tmp,char(176),Tbchar);
hold on
box on
axis square
set(gca,'fontsize', 15,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    %         if ischar(layerplot_color{layerplotnum})==1
    %             layerplot_color_here = char(layerplot_color(layerplotnum));
    %         else
    %             layerplot_color_here = layerplot_color{layerplotnum};
    %         end
    layerplot_color_here = 'k';
    layerplot_marker = '';
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)
    
    CPX=eval(sprintf('CS_minerals_OUT_%g(:,:,1)', layerplotnum));
    %CPX_1 = CS_minerals_OUT_1(:,:,1);
    % Ti is 4
    % Tb is 26
    plot(CPX(:,8), CPX(:,4)./CPX(:,8),'.k')
    plot(CPX(end,8), CPX(end,4)./CPX(end,8),'*r')
    
    %plot(eval(sprintf('CumFrac_%g', layerplotnum)),eval(sprintf('test_%g(:,10)', layerplotnum)),char(layerplot_marker(layerplotnum)),'Color',layerplot_color_here,'linewidth',3)
end
xlabel('Zr')
ylabel('Ti/Zr')
set(gca,'ticklength',2*[.02 .05])
%axis([.01 100 100 10000])
% set(gca,'yscale','log')
% set(gca,'xscale','log')
set(gcf,'Name','Zr_v_TiZr_in_CPX')



clear layerplotnum
figure('units','normalized','outerposition',[0.3935 0.3238 0.3923 0.5314])
%FigureTitle = sprintf('Ti Zr Tmp: %d %cC, %s', Tmp,char(176),Tbchar);
hold on
box on
axis square
set(gca,'fontsize', 15,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    %         if ischar(layerplot_color{layerplotnum})==1
    %             layerplot_color_here = char(layerplot_color(layerplotnum));
    %         else
    %             layerplot_color_here = layerplot_color{layerplotnum};
    %         end
    layerplot_color_here = 'k';
    layerplot_marker = '';
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    %         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)
    
    CPX=eval(sprintf('CS_minerals_OUT_%g(:,:,1)', layerplotnum));
    %CPX_1 = CS_minerals_OUT_1(:,:,1);
    % Ti is 4
    % Tb is 26
    plot(CPX(:,26), CPX(:,19)./CPX(:,26),'.k')
    plot(CPX(end,26), CPX(end,19)./CPX(end,26),'*r')
    
    %plot(eval(sprintf('CumFrac_%g', layerplotnum)),eval(sprintf('test_%g(:,10)', layerplotnum)),char(layerplot_marker(layerplotnum)),'Color',layerplot_color_here,'linewidth',3)
end
xlabel('Yb in CPX')
ylabel('Sm/Yb in CPX')
axis([1 20 10^-2 10])
set(gca,'ticklength',2*[.02 .05])
set(gca,'yscale','log')
set(gca,'xscale','log')
set(gcf,'Name','Yb_v_SmYb_in_CPX')


%%

clear layerplotnum
figure(); %('units','normalized','outerposition',[.5 .5 .8 .8])
FigureTitle = sprintf('Pressure v Bulk Compositions');
hold on

subaxis(1,3,1)
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    %         if ischar(layerplot_color{layerplotnum})==1
    %             layerplot_color_here = char(layerplot_color(layerplotnum));
    %         else
    %             layerplot_color_here = layerplot_color{layerplotnum};
    %         end
    layerplot_color_here = 'k';
    layerplot_marker = '';
    plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)
    
%     if layerplotnum==3
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color','g','MarkerSize',10)
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color','g','MarkerSize',10)
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color','g','MarkerSize',10)
%     end
    
    if layerplotnum==1
        leg=legend('First Melt Pressure','First Melt Extraction Pressure','Last Melt Pressure','autoupdate','off');
        set(leg,'Position',[0.1867 0.8744 0.1476 0.1083]);
%         
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color','b','MarkerSize',10)
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color','b','MarkerSize',10)
%         plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color','b','MarkerSize',10)
    end
    %plot(eval(sprintf('CumFrac_%g', layerplotnum)),eval(sprintf('test_%g(:,10)', layerplotnum)),char(layerplot_marker(layerplotnum)),'Color',layerplot_color_here,'linewidth',3)
    
    
end

axis([0.86 0.92 0 35])
axis square
% plot([0.02 0.02],[ 0 30],'r-')
%[Xtext,Ytext] = oxideLabel(1,Indicies(i));
xlabel('Mantle Mg# = MgO  / (MgO + FeO) in moles')
ylabel('Pressure (kbars)')
set(gca,'YDir','Reverse')
yline(Gar2SpTransition,'k');
yline(Sp2PlagTransition,'k');
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
%set(gca, 'Ytick',0:1:15)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.005:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):1:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
%%title(FigureTitle)


clear layerplotnum
% figure();%'units','normalized','outerposition',[.5 .5 .8 .8])
%
% FigureTitle = sprintf('Pressure Bulk NaK Tmp : %d %cC, %s', Tmp,char(176),Tbchar);
subaxis(1,3,2)
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
 
    layerplot_color_here = 'k';
    layerplot_marker = '*';
    plot(eval(sprintf('initialBulkComp_%g(12)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(12)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(12)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)

    
    
end

axis([0.04 0.20 0 35])
set(gca,'YDir','Reverse')
axis square
% plot([0.02 0.02],[ 0 30],'r-')
%[Xtext,Ytext] = oxideLabel(1,Indicies(i));
xlabel('Mantle NaK# = (Na_2O + K_2O) / (CaO + Na_2O + K_2O) in wt%')
ylabel('Pressure (kbars)')
yline(Gar2SpTransition,'k');
yline(Sp2PlagTransition,'k');
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
set(gca, 'Xtick',0.04:0.02:0.2)
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):1:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
%%title(FigureTitle)



subaxis(1,3,3)
axis square

clear layerplotnum
% figure();%'units','normalized','outerposition',[.5 .5 .8 .8])
% FigureTitle = sprintf('Pressure Bulk Ti Tmp : %d %cC, %s', Tmp,char(176),Tbchar);

hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
   
    layerplot_color_here = 'k';
    layerplot_marker = '*';
    
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),eval(sprintf('FirstMeltPb_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),eval(sprintf('FirstAt2percent_%g', layerplotnum)),'o','Color',fireEngineRed,'MarkerFaceColor',fireEngineRed,'MarkerSize',10)
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),eval(sprintf('LastMeltPb_%g', layerplotnum)),'o','Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',10)
    
    
    
end
axis([0 0.3 0 35])
set(gca,'YDir','Reverse')
% plot([0.02 0.02],[ 0 30],'r-')
%[Xtext,Ytext] = oxideLabel(1,Indicies(i));
xlabel('Mantle TiO_2 wt%')
ylabel('Pressure (kbars)')
yline(Gar2SpTransition,'k');
yline(Sp2PlagTransition,'k');
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
set(gca, 'Xtick',0:.05:.3)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):1:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
%%title(FigureTitle)

set(gcf,'position',[45 257 1233 420])



%%

clear layerplotnum
zeta = 100.*MASTER_ALL_RESULTS(:,19+4);
TMP =  MASTER_ALL_RESULTS(:,4);

figure(); 
FigureTitle = sprintf('Zeta v Bulk Compositions');
hold on

subaxis(1,3,1)
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
 
    layerplot_color_here = 'k';
    layerplot_marker = '';
    
    
    TMPhere = TMP(layerplotnum);
    P_index = find(TMPhere==TallOptions);
if isnan(zeta(layerplotnum))~=1
    plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),zeta(layerplotnum),'o','color',Tcolors{P_index},'MarkerSize',10)
end
    
   
    
end

axis([0.86 0.92 0 100])
axis square

xlabel('Mantle Mg# = MgO  / (MgO + FeO) in moles')
ylabel('\zeta in %')
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
%set(gca, 'Ytick',0:1:15)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.005:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)



clear layerplotnum
subaxis(1,3,2)
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    layerplot_color_here = 'k';
    layerplot_marker = '*';
 
    
    TMPhere = TMP(layerplotnum);
    P_index = find(TMPhere==TallOptions);
    if isnan(zeta(layerplotnum))~=1
    plot(eval(sprintf('initialBulkComp_%g(12)', layerplotnum)),zeta(layerplotnum),'o','color',Tcolors{P_index},'MarkerSize',10)
    end
    
end

axis([0.04 0.20 0 100])
axis square
% plot([0.02 0.02],[ 0 30],'r-')
%[Xtext,Ytext] = oxideLabel(1,Indicies(i));
xlabel('Mantle NaK# = (Na_2O + K_2O) / (CaO + Na_2O + K_2O) in wt%')
ylabel('\zeta in %')
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
set(gca, 'Xtick',0.04:0.02:0.2)
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
%%title(FigureTitle)



subaxis(1,3,3)
axis square
clear layerplotnum
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    
    layerplot_color_here = 'k';
    layerplot_marker = '*';
 
    TMPhere = TMP(layerplotnum);
    P_index = find(TMPhere==TallOptions);
    if isnan(zeta(layerplotnum))~=1
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),zeta(layerplotnum),'o','color',Tcolors{P_index},'MarkerSize',10)
    end
    
    
end
axis([0 0.3 0 100])
ylabel('\zeta in %')
xlabel('Mantle TiO_2 in wt%')
box on
grid on
%set(gca, 'Xtick',1500:10:1540)
set(gca, 'Xtick',0:.05:.3)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
%%title(FigureTitle)

set(gcf,'position',[45 257 1233 420])






%%
clear layerplotnum
figure()%'units','normalized','outerposition',[.5 .5 .8 .8])
FigureTitle = sprintf('CrustalThickness v Bulk Composition');
set(gcf,'position',[45 257 1233 420])

subaxis(1,3,1)
axis square
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
    
    layerplot_color_here = 'k';
    layerplot_marker = 'o';
    
    
    plot(eval(sprintf('initialBulkComp_%g(10)', layerplotnum)),eval(sprintf('CrustalThickness_%g(1)', layerplotnum)),char(layerplot_marker),'Color',layerplot_color_here,'MarkerSize',10,'MarkerFaceColor','k')
  
    
end

%set(gca,'YDir','Reverse')
axis([0.86 0.92 0 11])
set(gca, 'Ytick',0:11)
% plot([0.02 0.02],[ 0 30],'r-')
%[Xtext,Ytext] = oxideLabel(1,Indicies(i));
xlabel('Mantle Mg# =MgO/MgO+FeO) in moles')
ylabel('Crustal Thickness (km) ')
box on
grid on

ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.005:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 14,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)




subaxis(1,3,2)
axis square
clear layerplotnum
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
  
    layerplot_color_here = 'k';
    layerplot_marker = 'o';
    plot(eval(sprintf('initialBulkComp_%g(12)', layerplotnum)),eval(sprintf('CrustalThickness_%g(1)', layerplotnum)),char(layerplot_marker),'Color',layerplot_color_here,'MarkerSize',10,'MarkerFaceColor','k')
    
    
    
  
    
    
end
axis([0.04 0.2 0 11])
xlabel('Mantle NaK# =(Na_2O+K_2O)/(CaO+Na_2O+K_2O)in wt%')
ylabel('Crustal Thickness (km) ')
box on
grid on
set(gca, 'Ytick',0:1:11)

ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
set(gca, 'Xtick',0.04:0.02:0.2)
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 14,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')


subaxis(1,3,3)
axis square
clear layerplotnum
hold on
for layerplotnum = 1:size(MASTER_ALL_RESULTS,1)
   
    layerplot_color_here = 'k';
    layerplot_marker = 'o';
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),eval(sprintf('CrustalThickness_%g(1)', layerplotnum)),char(layerplot_marker),'Color',layerplot_color_here,'MarkerSize',10,'MarkerFaceColor','k')
    
    plot(eval(sprintf('initialBulkComp_%g(2)', layerplotnum)),eval(sprintf('CrustalThickness_%g(1)', layerplotnum)),char(layerplot_marker),'Color',layerplot_color_here,'MarkerSize',10,'MarkerFaceColor','k')
    
    
 
    
end
axis([0 0.3 0 11])
xlabel('Mantle TiO_2 wt%')
ylabel('Crustal Thickness (km) ')
box on
grid on
set(gca, 'Ytick',0:11)
set(gca, 'Xtick',0:.05:.3)

ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.5:ax.YLim(2);
set(gca,'ticklength',2*get(gca,'ticklength'))
set(gca,'fontsize', 14,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')


%%

figure();
FigureTitle = sprintf('Mean F versus Bulk Composition');
set(gcf,'position',[45 257 1233 420])
BulkNum = MASTER_POOLS_INFO(:,6);
BulkCompositionsNormalized = sumMgNumNormNoPT(BulkCompositions);
Comp2Plot = BulkCompositionsNormalized(BulkNum,:);


subaxis(1,3,1)
axis square
hold on
plot(Comp2Plot(:,11), MASTER_POOLS_INFO(:,29),'ko','MarkerSize',12)
axis([0.86 0.92 0 0.25])
xlabel('Mantle Mg#=MgO/(MgO+FeO) in moles')
ylabel('Fmean Full Pool')
box on
grid on
set(gca, 'Ytick',0:.05:.2)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.005:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.01:ax.YLim(2);
set(gca,'ticklength',2*[.02 .05])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)

subaxis(1,3,2)
axis square
hold on
axis([0.0 0.3 0 0.25])
plot(Comp2Plot(:,12), MASTER_POOLS_INFO(:,29),'ko','MarkerSize',12)
xlabel('Mantle NaK# = (Na_2O + K_2O) / (CaO + Na_2O + K_2O) in wt%')
ylabel('Fmean Full Pool')
box on
grid on
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.01:ax.YLim(2);
set(gca,'ticklength',2*[.02 .05])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)
set(gca, 'Xtick',0.0:0.05:1)
set(gca, 'Ytick',0.0:.05:.3)
%%title(FigureTitle)



subaxis(1,3,3)
axis square
hold on
%plot(MASTER_INSTANTS_INFO(:,5), MASTER_INSTANTS_DISEQ(:,1),'k.')
plot(Comp2Plot(:,2), MASTER_POOLS_INFO(:,29),'ko','MarkerSize',12)
axis([0 0.3 0.0 0.2])
xlabel('Mantle TiO_2 wt%')
ylabel('Fmean Full Pool')
box on
grid on
set(gca, 'Ytick',0:.05:.2)
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XLim(1):0.01:ax.XLim(2);
ax = gca;
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YLim(1):0.01:ax.YLim(2);
set(gca,'ticklength',2*[.02 .05])
set(gca,'fontsize', 13,'LineWidth',1,'FontName','Times New Roman')
set(gca,'XColor', 'k')
set(gca,'YColor', 'k')
set(gcf,'Name',FigureTitle)




